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In this paper we describe a method of calculation of master integrals 
based on the solution of systems of difference equations in one variable. 
Various explicit examples are given, as well as the generalization to arbi- 
trary diagrams. 

PACS numbers: ll.15.Bt; 02.90.+p 

1. Introduction 

Calculating Feynman integrals is not easy, especially when diagrams 
with many loops are considered. It is by now usual to use integration by 
parts [1, 2] identities for reducing a given generic Feynman integral to a 
linear combination of a small number of irreducible integrals, the so-called 
"master integrals". For these an explicit calculation is required. In this 
paper we illustrate a recently developed method of calculation based on the 
solution of difference equations. An extended exposition of this method 
and its applications can be found in Refs.[3, 4, 5, 6]. We mention here other 
different approaches with involve difference equations for calculating one- 
loop massive diagrams[7, 8], and structure functions for massless diagrams[9, 
10]. 



2. Difference equations 

Probably the reader is more accustomed to objects like differential equa- 
tions, which express relation between functions and their derivatives. Dif- 
ference equations are a particular kind of functional equations which express 
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relation between values of functions with argument shifted by integers. For- 
mally they are identical to "recurrence relations" but, instead of being only 
a mean to calculate the value of a function from its neighborhoods, they 
are functional equation in an unknown function; in some sense they may 
be considered as a "discretized" form of differential equations. In particu- 
lar, we will encounter linear difference equations. Difference equation are 
a well-known mathematical topic (G. Boole ~1850, N.E. Norlund 1929, 
L.M. Milne-Thomson 1933). A good exposition of the theory and methods 
of solution can be found in the beautiful book [11]. 

2.1. Fibonacci's numbers 

Let us start with a very basic example: the sequence of Fibonacci's inte- 
gers /(n) = {1, 1, 2, 3, 5, 8, 13, 21, . . .}; these numbers satisfy the recurrence 
relation 

/(n + 2) =/(n+l)+/(n) . (1) 

Now let us think of this relation as a seconG?-ordcr linear difference equation 
with constants coefficients. The general solution of this equation has the 
form I{n) = fJL^; by substituting it into the recurrence relation one obtains 
the characteristic equation 

/x^ - /X - 1 = , 

from which we work out the values of fi 

1 + V5 1-V5 

Therefore, the general solution of Eq.(l) is an arbitrary linear combination 
of the elementary solutions and order to reproduce the Fibonacci's 

sequence of integers we need to fix the constants, by using the initial con- 
ditions 7(1) = 7(2) = 1 . One finds 

3. One-loop vacuum integral 

Now let us jump directly to the calculation of the simplest master inte- 
gral: the euclidean integral in D dimensions (with unity mass for simplicity) 

n 
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(2) 

By using the identity obtained by integrating by parts [1, 2] 

[k-^ 1) 

one finds the recurrence relation 



(»-l) (^-(«-l-B/2) =0, (3) 

n n - 1 

which is a first-order difference equation with polynomial coefficients for the 
function J. A solution of Eq. (3) can be written in the form of an expansion 
in factorial series: 

h 'nn-K + s + 1) ^ 

— it^ r(n+l) , ai _| 02 1 

~ ^ r{n-K+l) "0 -r n-K+l ^ {n-K+l){n-K+2) ^ • " " 

where /i, K, Og are to be determined. 

3.0.1. Factorial series 

Factorial scries play for difference equations the same role that power 
series play for differential equations [11]. Factorial series of first kind like 
Eq.(4) are similar to asymptotic expansions in 1/n, but with the big advan- 
tage of be convergent! For example, let us consider the second derivative of 
the Gamma function 

*'(n)-^lnr(n). 

It satisfies the nonhomogeneous first order difference equation 

*'(n + l)-*'(n) = . 

If \I''(n) is expanded in power series of 1/n, one obtains an asymptotic series, 
not convergent for any value of n, 

where B2k are the Bernoulli's numbers; but if one expands it in factorial 
series, one obtains a series comfortably convergent for n > 0: 

y- r(^) r(n) 
^ 2^ s r(n + s) • 
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3.1. Solution of the difference equation 

By substituting the expansion (4) into the difference equation (3) (by 
using the Boole's operators[ll]) one finds the values = 1 and K = —D/2, 
and the recurrence relation between the coefficients 

sas = as-i{s + D/2){s + D/2 - 1) . 

The first coefficient ao is determined by comparing the large-n behaviours 
of the integral (2) and the factorial series (4): 



J (A;2 + l)n J 

asT{n + 1) 
T{n-K + s + l) 



s=0 

For s — > oo the term of the factorial series behaves as 

n el e^-l 

a.s S. S _ i+£,/2_„ . 

i> 1 



T{n + D/2 + s + l) .s! .s-D/2+n 
n = D/2 is the abscissa of convergence of the factorial series, in fact: 

• if n > D/2 the series is convergent, and it converges quickly forn 3> 1; 

• if n < Z?/2 the scries does not converge, but we can calculate J{n + i) 
for some large integer i by using the series, and then to use repeatedly 
the recurrence relation (3) in order to obtain J(n + i — J(n + i — 2), 

J(n). 

4. On-mass-shell self-mass master integral 

J {k^ + i)n(k^ -2p-k) " \y' 

n 

By combining identities obtained integrating by parts one finds a nonhomo- 
geneous second order difference equation 

{n-D)-Q)-+{2n-D-l)-Qy-3{n-l)-(y =(l-f) ^ . (5) 
n — 2 n — 1 n n — 1 

The general solution of this equation has the form 

7(n) = Ci7i(n) + C2/2(n) + hin) , 
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where Ii and I2 arc solutions of the homogeneous equation, Ci and C2 
arbitrary constants, and is a particular solution of the nonhomogeneous 
equation. As before, we look for solutions in the form of expansions in 
factorial series 

Un)-u-T ^-''^^(" + ^) ,--12 3 (6) 

s=Q ^ ■> ' 

By substituting Eq.(6) into Eq.(5) one finds the values 

/wi =1, M2 = -l/3, =1, 

= - D/2 + 1/2, K2=- D/2 + 1/2, K3 = - D/2 , 

and the recurrence relations 

4,fe(i) _ _ i)(7,- _D + + 3(.- - !)(.- - \fh% = , (7) 



(4s + 2)6f)_^(7s-L> + 4)6(-!^, + 3s-2(s- 1)622 = (f - l)(a,-sa,_i) , (8) 

where s = s + y — 1, aj = 6j = for f < and feQ"*^^ = h^^^ = 1. 

The constants Ci and C2 are determined by comparing the behaviours 
for n ^ 00, Ii(n) ^ n-^/2+V2^ /^(n) « (_i/3)n„-D/2+i/2 ^^^^ j^^^) ^ 

(l/2-D/4)n--°/2 with 

V«y J k^-2p-k 2 

n 

one finds Ci = C2 = 0, therefore the solution I2 can be discarded. As 
in the previous case, the factorial series converge quickly for ra 3> 1, and 
the recurrence relations are used for the values of n such that the series 
converges slowly or does not converge. 

5. General case 

Now we sketch the generalization of the method used above to the calcu- 
lation of master integrals of arbitrary diagrams. Given an L-loop integral, 
we put the exponent n on the first denominator (arbitrarily chosen) 
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the numerator N in general contain products of powers of some irreducible 
scalar products involving loop momenta. By combining by suitable algo- 
rithms the identities obtained by integration by parts 



/ 



one finds a nonhomogeneous difference equation of order r (typically r ~ 
2-10) 

Po{n)B{n)+pi{n)B{n + l) + ...+Pr{n)B{n + r) = T{n) ; 

the right-hand side T(n) contains new integrals Bi, B2,. . . which have one 
or more denominators missing with respect to as an example 

N 



'D^D^D^...Dn, 

and which are simpler. By building difference equations for ^2, . . . one 
obtains a system of difference equations in triangular form. The system 
is solved one equation at once, beginning with that corresponding to the 
simplest master integral, with the minimum number of denominators. The 

solutions 

Bi{n) = Ci/f ^(n) + C2/2°^(n) + . . .^/.^^(n) + /^^(n) 
are expanded in series 



tOM( . _ n sr^ r(n + l) 

^3 A.r(n-i^,+s + l) 



s= 

S 00 



fJ'f'^, and ai^'^'^^ are obtained from the expansion of T{n) which 

is supposed to be known. By substituting the expansions into the dif- 
ference equation one finds the characteristic equation and the values of 
fXi, fi2, ■ ■ ■, l-i-r, the r indicial equations and the values of Ki,. . . , Kr, and the 
recurrence relations for the coefficients which have usually order higher 
than Eqs.(7)-(8)(~10-20). By comparing the behaviours for n ^ 00 the 
values of the constants Ci, C2, ■ ■ ■ , Cr arc obtained (usually most of con- 
stants are null). For example, the large-n behaviour of an euclidean massive 
integral is: 

/ N n-D/2 
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where Di = kf + mf, mi ^ 0. F{0) is the value of the integral corresponding 
to the diagram with one loop less obtained by setting A;i=0, graphically 




n 



6. System of difference equations in the single-scale case 

To give to the reader an idea of how the difference equations for more 
complicated diagrams look, we list here the whole system of difference equa- 
tions for the two-loop self-mass diagram in the particular case with all 
masses equal to 1 and on-shell external momentum. As above, the point 
on a line means that the corresponding denominator is raised to the power 
indicated. The factorizable two-loop subdiagrams are decomposed into the 
product of one-loop diagrams. 

^ _(„_i_£,/2) ^ =0, 
n n-1 



(„_Z))_Q_+(2n-i?-l)-Q--3(n-l)-Q- = (1-f ) ^ , 



n- 2 



n-1 



n-1 



^n-D) Q +{2n-D-l) Q -3(n-l) Q = {2-D) x () , 



n-2 



n-1 



n-1 



{n-D){2n-?,D + 2) +2[7r? - n{lQD + 1) + 2,0'^ + -2] 

n — 2 n — 1 

-16(n-l)(n-D+l) = 2{D-2f x , 



n-1 



n-2 



n-1 



(2n-3L>+4) -Q- +(4n-D-2) -Q- -6(n-l) -Q- = 
-(D-2)Q -2(D-2)-Q- X 



n-1 



n-1 
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48(n-D + l) 



48(n-D + l) 



-i^)-^+(2n-i5-l)-Q--3(n-l)-^j = 
n — 2 n — 1 n 

-5(n-r>)(2n-3L> + 2) -Q- - [-lOn^ + n(28D - 58) - 18D^ 

n-2 

+ 46D + 12]-Q- -(£>-2)(-32n + 6D + 52) x ^ . 
n — 1 n — 1 

n — 2 n — 1 n -. 

(^n-D)-^y +{2n-D-l)-^y -3(n-l)-^^ = 

-5(n-i:')(2ra-3i:' + 2) -[-lOn'^ + n{28D-58)-18D'^ + 

n-2 

+ 12] -Q- -24(L>-2)(n-L> + l) ^ x -Q- 
n- 1 

-2(D-2)(8n-9D+2) x ^ 



n - 1 



n - 1 
n-2 



n- 1 



36(n-L> + l) 2(n-D + l) -(J)- + (2n-D) -^J)--4(n-l) -(J)- = 

[n — 2 n — 1 

(2n - 3£> + 4) -Q- - 2(n - D + 1) -Q- 

n - 2 

+ (n-2£> + 4) -(n-i:>)(2ri-3£) + 2) -Q- 

n - 2 

- (-2n2 - 16ni:> + 46n + ISi:*^ - 70D + 60) + 12(n - D + 1) |^ 

n — 1 

-(n-L')-Q- X -Q- +(n-2L> + 4) -Q- x -Q- 

n — 2 n — 1 

+ (0-2) -p-2) (3 X -Q-(0-2)-(;> X ^] 

n — 1 n — 1 

+ 2(L>-2)(8n-3L>-10) ^ ^ . 



n- 1 



n-1 
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7. The program SYS 

A generic integral, through the solution of a system of algebraic identities 
between integrals obtained integrating by parts, is transformed into a linear 
combination of a small number of "master integrals" , chosen typically with 
numerator unity and denominators raised to power one. 

The method of calcTiIation based on the solution of finite difference equa- 
tions is applied only to these "master integrals" . 

Considered the number and the complexity of the systems of difference 
equations for diagrams with more than 2 loops, we have written on purpose 
a comprehensive program, called SYS, with the following features: 

• C program, about 23000 lines. 

• The program automatically determines the master integrals of a dia- 
gram, it builds and solves the systems of difference equations. 

• Input: description of the diagram, number of terms of the expansion 
in D - 4. 

• The program contains a simplified algebraic manipulator, used to solve 
systems of identities among integrals with this kind of coefficients: 
arbitrary precision integers, rationals, ratios of polynomials in one 
and two variables (for example D and x) with integer coefficients. 

• Efficient management of systems of identities of size up to the limit of 
disk space (tested up to 5 • 10^ identities). 

• Numerical solution of systems of difference and differential equations 
up to 500 equations, using arbitrary precision floating point complex 
numbers and truncated series in e. 

• All the coefficients of the expansions in e are worked out in numerical 
form, even those of divergent terms. 

• Floating number precision: up to 1200 digits (see [6]) (essentially one 

sums expansions in one variable). 

• Arithmetic libraries which deal with operations on integers, polyno- 
mials, rationals, floating point numbers and truncated series in e were 
written on purpose by the author. 

8. Applications of the method 

1. high-precision calculation of vacuum (2-4 loops), self-mass (1-3 loops), 
vertex and box diagrams (1-2 loops) in the case of equal masses and 
on-shell external momenta [3, 4, 6]. 
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2. high-precision calculation of electron g — 2 master integrals, at 3 loops 
[5] and at 4 loops[12] (actually only a few, work in progress). 

3. high-precision calculation of integrals for arbitrary values of mass and 
momenta (work in progress). 

Here is an example of high-precision (60 digits) result obtained with the 
program SYS: the expansion in e = (4 — D)/2 of the following two-loop self- 
mass integral with these values of masses and momenta: = —1, ml = 2, 
m| = 3, m| = 4, ml = 5, m^ = 6, 




= r(l + e)^[ 



+ 0.213887190522735619987564582191165767549812066097826018573 
-0.918487358981510992237303382863013109752582619839173727799e 
-|-2.395114573687568040190458275515297803861359807588836500079e^ 
- 5.204471820133273025267916826530760950863731824065596549403e^ 
-M0.62755429637642806105390035038189357736980511197603387938e*-hO(e^)] . 



9. Conclusions 

1. Method of calculation applicable to diagrams with any topology. 

2. The calculation of an integral is reduced to the sum of convergent series 
in one variable. As consequence, the number of operations needed to 
obtain a precision of the results of N digits grows linearly with N and 
does not depends on the complexity of the diagram; calculations with 
N ~ 10^ are still feasible[6]. 

3. Easiness of cross-checks, due to possibility of writing different differ- 
ence equations for each topologically distinct line of a given diagram. 
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